% This script starts estimation of the Macro finance model using the SR approach
% By Martin M. Andreasen, August 2019
addpath(genpath(pwd))
rmpath([pwd,'\SRstep2_threshold'])
close all
clear 
%clc
%% User settings
dataOption        = 1;          %1 for using ygap = UGAP,
                                %2 for using ygap = potential output in GDP
                                %3 for using ygap = Cooper and Priestley gap
startYear         = 1961;       %startYear for the estimation. Data available from 1961
modelNum          = 1;          %1 for QTSM without interest rate smoothing
                                %2 for QTSM wiht constant interest rate smoothing

% For inference
theta12OptimalOn  = 0;          %1 for setting LAMBDA in step 3 in an optimal manner. 0 for using the estimates from step 2
AVarTheta1On      = 0;          %1 for computing AVar for theta1 at step 1 and theta11 at step 3, else 0

% The optimizer:
optimizerStep1     = 3*0;         % optim = 0 for no optimization
optimizerStep3     = 3*0;         % optim = 1 for the CMAES
                                % optim = 2 for a gradient based optimizer (the LM version)
                                % optim = 3 for the Nelder-Mead.
                                % optim = 4, for first a gradient optimizer and the the Nelder-Mead
optim.numOptim    = 4;          % Number of optimizations - if optimizer = 3 or 4.
optim.sigma       = 0.10;       % step length for CMAES
optim.sigmaMax    = 0.20;       % max length for CMAES
optim.TolFun      = 1e-5;
optim.TolX        = 1e-5;
optim.MaxIter     = 4000;
optim.MaxFunEvals = 4000;

epsValue          = 1D-5;       % Stepsize for the standard errors of Theta1
AVarAdjThetaOn    = 0;          % 1 to adjust factor variance for Theta1, else 0

h0RegimeOn        = 0;          %1 for switching in h0 under P, else 0
hxRegimeOn        = 1;          %1 for switching in hx under P, else 0

%% Default setting
numCPUs           = feature('numcores'); % Number of CPUs
fullPHIon         = 1;          %1 for a full phi, else 0 for diagonal phi
nm                = 2;          %Number of macro factors
nn                = 0;          %Number of latent factors factors
loadDefaultSetting
%factorSignResOn = 2

% Loading the data:
loadDataAll = getData(dataOption,startYear);
% Reduce sample to end at the end of 2007. 
datestr(loadDataAll.yieldsOut(2:560,1));
loadData.yieldsOut       = loadDataAll.yieldsOut(1:561,:);
loadData.yields_h1       = loadDataAll.yields_h1(1:561,:);
loadData.crspRiskFree_h1 = loadDataAll.crspRiskFree_h1(1:560,:);
loadData.crspRiskFree_h3 = loadDataAll.crspRiskFree_h3(1:560,:);
loadData.NBERrec         = loadDataAll.NBERrec(1:560,:);
loadData.infl            = loadDataAll.infl(1:560,:);
loadData.ygap            = loadDataAll.ygap(1:560,:);
loadData.econ            = loadDataAll.econ(1:560,:);
%% Starting values for the estimated parameters
theta1 = initializeTheta1(fullPHIon,nm,nn);

% Impose restrictions on phi from ReStud paper
theta1 = rmfield(theta1,'phi14');
theta1 = rmfield(theta1,'phi23');
theta1 = rmfield(theta1,'phi32');
theta1 = rmfield(theta1,'phi34');
theta1 = rmfield(theta1,'phi41');
theta1 = rmfield(theta1,'phi43');

if modelNum == 1  %
    theta1 = rmfield(theta1,'rhor');
    theta1ValuesStep1 = [  0.000002650478    0.089197446908    0.000001021766    0.006343174840    0.025203528490  ...
        0.003735193053    0.003955627181   -0.488101502565   -0.488047311878   -0.023233075239  ...
        -0.000900122006    0.002768797685   -0.023763822785    0.000225891899    0.000837232199  ...
        0.023921163031   -0.089100997670    0.049999963490    0.018081454574    0.000990441851  ...
        0.041104642628    0.002548114235   -0.004158987706    0.002122385727    0.000002404915  ]';
    % Q = 0.067669314609665, 
    
    % For step 3
    if h0RegimeOn == 0 && hxRegimeOn == 0
        theta1ValuesStep3 = [    0.000151995377    0.062085877129    0.000001393510    0.023306719513    0.022570605271  ...
            0.002911844185    0.003612797338   -0.252180756766   -0.333831082668    0.025039585434  ...
            -0.000925302090    0.002189387848   -0.014186606604    0.005644979139    0.000554509835  ...
            0.001742256447    0.000061441755    0.002035395934   -0.006870024697   -0.002088868750  ...
            0.012660665036   -0.000769418668   -0.000215944113    0.000978490257    0.000690464334  ]';
        % Q =   0.070302926395732, done
    elseif h0RegimeOn == 1 && hxRegimeOn == 0
        theta1ValuesStep3 = [ ]';
        % Q =
    elseif h0RegimeOn == 0 && hxRegimeOn == 1
        theta1ValuesStep3 = [   0.000180657136    0.061942212830    0.000001687429    0.023153885926    0.022243584067  ...
            0.002836286665    0.003446942321   -0.250433238530   -0.322512789508    0.026191460396  ...
            -0.000919395196    0.002155789239   -0.013992850656    0.005654175469    0.000550819707  ...
            0.001781401355    0.000065916619    0.002033390898   -0.006746296841   -0.001948361634  ...
            0.012539790957   -0.000728960237   -0.000209971531    0.001007439075    0.000701314765  ]';
        % QStep3 = 0.070271994613864, done
    elseif h0RegimeOn == 1 && hxRegimeOn == 1
        
    end
elseif modelNum == 2
    
end
theta1           = values2struct(theta1ValuesStep1,fieldnames(theta1));
theta1StartStep3 = values2struct(theta1ValuesStep3,fieldnames(theta1));

% Getting bounds and Insigma
[upperBounds,lowerBounds,Insigma] = theta1BoundsInsigma(nm,nn);

% Defining the struct with the setup specification for step 1
constructSetupStep1

%% Computing the objective function theta1
tic
[Q,output,errorMes] = objectFunc(struc2values(theta1),setupStep1);
toc
%% Step 1 of the SR approach
% The optimization
optim.optimizer = optimizerStep1;
[QStep1,theta1Step1,outputStep1] = step1SR(theta1,setupStep1,optim);
theta1 = theta1Step1;

% The standard errors
AVarTheta1Step1  = getAVarTheta1andFactors(theta1Step1,setupStep1,epsValue);
%plotFactors(setupStep1,outputStep1,AVarTheta1Step1)
%% Step 2 of the SR approach
[theta2Step2,AVarTheta2Step2] = step2SR_threshold_macro(AVarTheta1Step1,outputStep1,bandwidthTstep2,setupStep1.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
%mprStep2 = getMPR(outputStep1,theta2Step2,AVarTheta2Step2,nm*2+nn,h0RegimeOn,hxRegimeOn);
%% Step 3 of the SR approach
namesTheta12 = cell((2*nm+nn)*(2*nm+nn+1)/2,1);
idx = 0;
for i=1:2*nm+nn
    for j=1:i
        idx = idx + 1;
        namesTheta12{idx} = ['sigma',num2str(i),num2str(j)];
    end
end
[theta12Step3,AVarTheta12Step3] = step3SR_theta12(namesTheta12,theta1Step1,AVarTheta1Step1,theta2Step2,AVarTheta2Step2,theta12OptimalOn);

% Estimation of theta11 - using Robust version
if exist('theta1StartStep3','var')
    theta1Start = theta1StartStep3;
    if AVarTheta1On  == 1
        disp('for standard errors')
        theta1Start = rmfield(theta1Start,'phi22');
        setupStep1.calibrateTheta1.phi22 = theta1StartStep3.phi22;
        setupStep1.selectTheta1          = fieldnames(theta1Start);       
        setupStep1.upperBoundsValues    = struc2values(upperBounds,setupStep1.selectTheta1);
        setupStep1.lowerBoundsValues    = struc2values(lowerBounds,setupStep1.selectTheta1);
        setupStep1.InsigmaValues        = struc2values(Insigma,setupStep1.selectTheta1);
     end
else
    theta1Start = theta1Step1;
end
optim.optimizer = optimizerStep3;
[QStep3,theta11Step3,outputStep3,setupStep3] = step3SR_theta11(theta12Step3.Robust,namesTheta12,theta1Start,...
    setupStep1,upperBounds,lowerBounds,Insigma,optim);

% Standard errors of theta11: The asymptotic standard errors 
AVarTheta11Step3 = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step2,AVarTheta12Step3,...
    setupStep1,setupStep3,epsValue/10,theta12OptimalOn);

% Re-estimating P dynamics
[theta2Step3,AVarTheta2Step3] =step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);

% Enforce consistency of sigma between Q and P
posTheta12Step2 = zeros(1,length(namesTheta12));
for i=1:length(namesTheta12)
    theta2Step3.(namesTheta12{i}) = theta2Step2.(namesTheta12{i});
    posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step3.names,namesTheta12(i)));
end
AVarTheta2Step3.VarRobust(posTheta12Step2,posTheta12Step2) = AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2); 
% Check the two estimates of sigma
[h0_1,h0_2,hx_1,hx_2,sigmaStep3] = unfoldTheta2_threshold_macro(theta2Step3,nm*2+nn);
%[h0_1,h0_2,hx_1,hx_2,sigmaStep2] = unfoldTheta2_threshold_macro(theta2Step2,nm*2+nn);

mprStep3 = getMPR(outputStep3,theta2Step3,AVarTheta2Step3,nm*2+nn,h0RegimeOn,hxRegimeOn);
testRegimeStep3 = waldTestRegimeShift(theta2Step3,AVarTheta2Step3,2*nm+nn);

%% Do the yield curve decomposition
if h0RegimeOn == 0 && hxRegimeOn == 0
    yieldTP = yieldCurveDecom(setupStep3,outputStep3,theta2Step3);
else
    numSimTP = 1D4;
    yieldTP  = yieldCurveDecom_threshold(setupStep3,outputStep3,theta2Step3,numSimTP);
end

% mean(yieldTP.termPremia(end,:))
% std(yieldTP.termPremia(end,:))
%% Simulating the model
thresholdLevel          = AVarTheta2Step3.thresholdLevel;
momentsData             = getMomentsData(loadData,outputStep3,thresholdLevel);
momentsModel            = simulator_threshold(theta2Step3,outputStep3,thresholdLevel,numSim);
momentsModel_finite     = simulator_threshold_finite(theta2Step3,outputStep3,thresholdLevel,1000);
expectedExcessReturn.h1 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,1);
expectedExcessReturn.h3 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,3);

%% Do plots
plotFactors(setupStep3,outputStep3,AVarTheta11Step3)
plotFit(setupStep3,outputStep3)
plotUnconditionalMoments(momentsModel,momentsData)
plotCSregressions(momentsModel,momentsData)
plotForecastRegressions(momentsModel,momentsData)
plotYieldCurveDecom(yieldTP,setupStep3)


%% Displaying the results
outTheta1          = displayResults(theta1Step1,AVarTheta1Step1.VarRobust);
outTheta11         = displayResults(theta11Step3,AVarTheta11Step3.VarRobust);
theta2Step3Reduced = reduceTheta2(theta2Step3,AVarTheta2Step3.names);
outTheta2Reduced   = displayResults(theta2Step3Reduced,AVarTheta2Step3.VarRobust);

disp('optimum at step 1')
printToScreen(struct2array(theta1Step1));
disp('optimum at step 3')
theta1Step3Start  = theta11Step3;
for i = 1:size(namesTheta12,1)
    theta1Step3Start.(namesTheta12{i}) = setupStep3.calibrateTheta1.(namesTheta12{i});
end
printToScreen(struct2array(theta1Step3Start));

%% Saving the results in a mat file
if modelNum == 1
    modelName = 'QTSM_nn0_ReStudPhiQ_end2007_';
elseif modelNum == 2
    modelName = 'QTSMr_nn0_ReStudPhiQ_end2007_';
end
filename = [modelName,num2str(startYear),'_SE',num2str(epsValue),'_h0RegimeOn',num2str(h0RegimeOn),'_','hxRegimeOn',num2str(hxRegimeOn),...
    '_numBootStep2_',num2str(numBootStep2)];
save([pwd,'\Results\',filename,'.mat'],'QStep1','theta1Step1','outputStep1','AVarTheta1Step1',...
    'theta2Step2','AVarTheta2Step2',...
    'theta11Step3','outputStep3','setupStep3','AVarTheta11Step3','theta2Step3','AVarTheta2Step3',...
    'yieldTP','momentsModel','momentsModel_finite','momentsData','mprStep3','testRegimeStep3','expectedExcessReturn')

